
logAdd <- function(lx, ly) {
  max(lx, ly) + log1p(exp(-abs(lx - ly)))
}


mstepMisreport <- function(y, x.misreport, w, treat,
                              include.misreport.treatment, weight) {

  lrep <- rep(c(1, 0), each = length(y)) 

  if(include.misreport.treatment == TRUE) {
    xrep <- as.matrix(rbind(cbind(x.misreport, treat), cbind(x.misreport, treat)))
  } else if(include.misreport.treatment == FALSE) {
    xrep <- as.matrix(rbind(x.misreport, x.misreport))
  }

  wrep <- c(w[, 1] + log(weight), w[, 2] + log(weight)) 

  lrep <- lrep[wrep > -Inf]
  xrep <- xrep[wrep > -Inf, , drop = FALSE]
  wrep <- wrep[wrep > -Inf]

  X <- xrep

  if(ncol(X) == 1) {
    fit.misreport <- glm(cbind(lrep, 1 - lrep) ~ 1,     weights = exp(wrep), family = binomial)
  } else if(ncol(X) > 1) {
    fit.misreport <- glm(cbind(lrep, 1 - lrep) ~ -1 + X, weights = exp(wrep), family = binomial)
  }

  coefs <- coef(fit.misreport)
  names(coefs) <- gsub("^X1|^X2|^X3|^X", "", names(coefs))

  return(coefs)

}


mstepSensitive <- function(y, treat, x.sensitive, w, d, sensitive.response,
                           weight, model.misreport) {

  if(model.misreport == TRUE) {
    zrep <- rep(c(sensitive.response, abs(1 - sensitive.response)), each = length(y))
    xrep <- as.matrix(rbind(x.sensitive, x.sensitive))
    wrep <- c(apply(w[, 1:2], 1, function(x) logAdd(x[1], x[2])) + log(weight), w[, 3] + log(weight))

    

    zrep <- zrep[wrep > -Inf]
    xrep <- xrep[wrep > -Inf, , drop = FALSE]
    wrep <- wrep[wrep > -Inf]

    X <- xrep

    if(ncol(X) == 1) fit <- glm(cbind(zrep, 1 - zrep) ~ 1,     weights = exp(wrep), family = binomial)
    if(ncol(X) > 1)  fit <- glm(cbind(zrep, 1 - zrep) ~ -1 + X, weights = exp(wrep), family = binomial)

    coefs <- coef(fit)
  } else {
    zrep <- rep(c(1, 0), each = length(y))
    xrep <- as.matrix(rbind(x.sensitive, x.sensitive))
    wrep <- c(w[, 1] + log(weight), w[, 2] + log(weight))

    zrep <- zrep[wrep > -Inf]
    xrep <- xrep[wrep > -Inf, , drop = FALSE]
    wrep <- wrep[wrep > -Inf]

    X <- xrep

    if(ncol(X) == 1) fit <- glm(cbind(zrep, 1 - zrep) ~ 1,     weights = exp(wrep), family = binomial)
    if(ncol(X) > 1)  fit <- glm(cbind(zrep, 1 - zrep) ~ -1 + X, weights = exp(wrep), family = binomial)

    coefs <- coef(fit)
  }

  names(coefs) <- gsub("^X1|^X2|^X3|^X", "", names(coefs))

  return(coefs)

}


mstepControl <- function(y, treat, J, x.control, w, d, sensitive.response,
                         weight, model.misreport, control.constrained) {

  if(model.misreport == TRUE) {

    if(control.constrained == "none") {
      yrep <- c((y - treat * as.numeric(sensitive.response == 1)),
                (y - treat * as.numeric(sensitive.response == 1)),
                (y - treat * as.numeric(sensitive.response == 0)))
      xrep <- as.matrix(rbind(x.control, x.control, x.control))
      zrep1 <- rep(c(1, 0, 0), each = length(y)) 
      zrep2 <- rep(c(1, 1, 0), each = length(y)) 
      zrep3 <- rep(c(0, 0, 1), each = length(y)) 
      wrep <- c(w[, 1] + log(weight), w[, 2] + log(weight), w[, 3] + log(weight))

      yrep <- yrep[wrep > -Inf]
      xrep <- xrep[wrep > -Inf, , drop = FALSE]
      zrep1 <- zrep1[wrep > -Inf]
      zrep2 <- zrep2[wrep > -Inf]
      zrep3 <- zrep3[wrep > -Inf]
      wrep <- wrep[wrep > -Inf]

      X <- cbind(xrep, U = zrep1, Z = zrep2)
      control.fit <- glm(cbind(yrep, J - yrep) ~ -1 + X, weights = exp(wrep), family = binomial)
      coefs <- coef(control.fit)

    } else if(control.constrained == "partial") { 
      yrep <- c((y - treat * as.numeric(sensitive.response == 1)),
                (y - treat * as.numeric(sensitive.response == 0)))
      xrep <- as.matrix(rbind(x.control, x.control))
      zrep1 <- rep(c(1, 0), each = length(y)) 
      zrep2 <- rep(c(0, 1), each = length(y)) 
      wrep <- c(apply(w[, 1:2], 1, function(x) logAdd(x[1], x[2])) + log(weight), w[, 3] + log(weight))

      yrep <- yrep[wrep > -Inf]
      xrep <- xrep[wrep > -Inf, , drop = FALSE]
      zrep1 <- zrep1[wrep > -Inf]
      zrep2 <- zrep2[wrep > -Inf]
      wrep <- wrep[wrep > -Inf]

      X <- cbind(xrep, Z = zrep1)
      control.fit <- glm(cbind(yrep, J - yrep) ~ -1 + X, weights = exp(wrep), family = binomial)
      coefs <- coef(control.fit)

    } else if(control.constrained == "full") { 
      yrep <- c((y - treat * as.numeric(sensitive.response == 1)),
                (y - treat * as.numeric(sensitive.response == 0)))
      xrep <- as.matrix(rbind(x.control, x.control))
      wrep <- c(apply(w[, 1:2], 1, function(x) logAdd(x[1], x[2])) + log(weight), w[, 3] + log(weight))

      yrep <- yrep[wrep > -Inf]
      xrep <- xrep[wrep > -Inf, , drop = FALSE]
      wrep <- wrep[wrep > -Inf]

      X <- xrep
      
      if(ncol(X) == 1) control.fit <- glm(cbind(yrep, J - yrep) ~ 1    , weights = exp(wrep), family = binomial)
      if(ncol(X) > 1)  control.fit <- glm(cbind(yrep, J - yrep) ~ -1 + X, weights = exp(wrep), family = binomial)
      coefs <- coef(control.fit)
    }
  } else {

      yrep <- c(y - treat, y)
      xrep <- as.matrix(rbind(x.control, x.control))
      zrep1 <- rep(c(1, 0), each = length(y))
      zrep2 <- rep(c(0, 1), each = length(y))
      wrep <- c(w[, 1] + log(weight), w[, 2] + log(weight))

      yrep <- yrep[wrep > -Inf]
      xrep <- xrep[wrep > -Inf, , drop = FALSE]
      zrep1 <- zrep1[wrep > -Inf]
      zrep2 <- zrep2[wrep > -Inf]
      wrep <- wrep[wrep > -Inf]

    if(control.constrained == "none") {
      X <- cbind(xrep, Z = zrep1)
      fit.partial <- glm(cbind(yrep, J - yrep) ~ -1 + X, weights = exp(wrep), family = binomial)
      coefs <- c(coef(fit.partial))
    }

    if(control.constrained == "full") {
      X <- xrep
      if(ncol(X) == 1) fit.full <- glm(cbind(yrep, J - yrep) ~ 1    , weights = exp(wrep), family = binomial)
      if(ncol(X) > 1)  fit.full <- glm(cbind(yrep, J - yrep) ~ -1 + X, weights = exp(wrep), family = binomial)
      coefs <- c(coef(fit.full))
    }
  }

  names(coefs) <- gsub("^X1|^X2|^X3|^X", "", names(coefs))
  if(model.misreport == TRUE && sensitive.response == 1) names(coefs)[names(coefs) %in% c("Z", "")] <- "Z == 1"
  if(model.misreport == TRUE && sensitive.response == 0) names(coefs)[names(coefs) %in% c("Z", "")] <- "Z == 0"
  if(model.misreport == FALSE) names(coefs)[names(coefs) %in% c("Z", "")] <- "Z == 1"
  names(coefs)[names(coefs) == "(Intercept):1"] <- "(Intercept)"

  return(coefs)

}


mstepOutcome <- function(y, treat, x.outcome, w, d, sensitive.response, o,
                         trials, weight, outcome.model, model.misreport,
                         outcome.constrained, control.constrained) {

  coefs.aux <- NULL

  if(outcome.constrained == TRUE) {
    if(model.misreport == TRUE) {
      xrep <- as.matrix(rbind(x.outcome, x.outcome))
      zrep <- rep(c(1, 0), each = length(y))
      orep <- as.matrix(c(o, o))
      trialsrep <- as.matrix(c(trials, trials))
      wrep <- c(apply(w[, 1:2], 1, function(x) logAdd(x[1], x[2])) + log(weight), w[, 3] + log(weight))
    } else {
      xrep <- as.matrix(rbind(x.outcome, x.outcome))
      zrep <- rep(c(1, 0), each = length(y))
      orep <- as.matrix(c(o, o))
      trialsrep <- as.matrix(c(trials, trials))
      wrep <- c(w[, 1] + log(weight), w[, 2] + log(weight))
    }

    xrep <- xrep[wrep > -Inf, , drop = FALSE]
    zrep <- zrep[wrep > -Inf]
    orep <- orep[wrep > -Inf]
    trialsrep <- trialsrep[wrep > -Inf]
    wrep <- wrep[wrep > -Inf]

    X <- cbind(xrep, Z = zrep)[, -1, drop = FALSE]

    if(outcome.model == "logistic") {
      fit.constrained.logistic <- glm(cbind(orep, 1 - orep) ~ 1 + X, weights = exp(wrep), family = binomial)
      coefs <- coef(fit.constrained.logistic)
    } else if(outcome.model == "binomial") {
      fit.constrained.binomial <- glm(cbind(orep, trialsrep - orep) ~ 1 + X, family = binomial, weights = exp(wrep))
      coefs <- coef(fit.constrained.binomial)
    } else if(outcome.model == "betabinomial") {
      fit.constrained.betabinomial <- VGAM::vglm(cbind(orep, trialsrep - orep) ~ 1 + X, VGAM::betabinomial, weights = exp(wrep))
      coefs <- coef(fit.constrained.betabinomial)[-2]
      coefs.aux <- c(rho = mean(fit.constrained.betabinomial@misc$rho))
    }
  } else if(outcome.constrained == FALSE) {
    if(model.misreport == TRUE) {
      xrep <- as.matrix(rbind(x.outcome, x.outcome, x.outcome))
      zrep1 <- rep(c(1, 0, 0), each = length(y))
      zrep2 <- rep(c(1, 1, 0), each = length(y))
      orep <- as.matrix(c(o, o, o))
      trialsrep <- as.matrix(c(trials, trials, trials))
      wrep <- c(w[, 1] + log(weight), w[, 2] + log(weight), w[, 3] + log(weight))
    } else {
      stop("\noutcome.constrained = TRUE is only possible when a direct question is included.")
    }

    xrep <- xrep[wrep > -Inf, , drop = FALSE]
    zrep1 <- zrep1[wrep > -Inf]
    zrep2 <- zrep2[wrep > -Inf]
    orep <- orep[wrep > -Inf]
    trialsrep <- trials[wrep > -Inf]
    wrep <- wrep[wrep > -Inf]

    X <- cbind(xrep, U = zrep1, Z = zrep2)[, -1, drop = FALSE]

    if(outcome.model == "logistic") {
      fit.unconstrained.logitistic <- glm(cbind(orep, 1 - orep) ~ 1 + X, weights = log(wrep), family = binomial)
      coefs <- coef(fit.unconstrained.logitistic)
    } else if(outcome.model == "binomial") {
      fit.unconstrained.binomial <- glm(cbind(orep, trialsrep - orep) ~ 1 + X, family = binomial, weights = log(wrep))
      coefs <- coef(fit.unconstrained.binomial)
    } else if(outcome.model == "betabinomial") {
      fit.constrained.betabinomial <- VGAM::vglm(cbind(orep, trialsrep - orep) ~ 1 + X, VGAM::betabinomial, weights = log(wrep))
      coefs <- coef(fit.constrained.betabinomial)[-2]
      coefs.aux <- c(rho = mean(fit.constrained.betabinomial@misc$rho))
    }
  }

  names(coefs) <- gsub("^X", "", names(coefs))
  names(coefs)[names(coefs) == ""] <- "Z"
  names(coefs)[names(coefs) == "(Intercept):1"] <- "(Intercept)"

  return(list(coefs = coefs, coefs.aux = coefs.aux))

}


estep <- function(y, w, x.control, x.sensitive, x.outcome, x.misreport, treat, J,
                  par.sensitive, par.control, par.outcome,
                  par.outcome.aux, par.misreport,
                  d, sensitive.response, model.misreport,
                  o, trials, outcome.model, weight,
                  outcome.constrained, control.constrained, respondentType,
                  include.misreport.treatment) {

  log.lik <- rep(as.numeric(NA), length(y))

  if(model.misreport == TRUE) {

    
    if(control.constrained == "none") {
      hX.misreport.sensitive <- plogis(cbind(x.control, 1, 1) %*% par.control, log.p = TRUE)
      hX.truthful.sensitive <- plogis(cbind(x.control, 0, 1) %*% par.control, log.p = TRUE)
      hX.truthful.nonsensitive <- plogis(cbind(x.control, 0, 0) %*% par.control, log.p = TRUE)
    }

    if(control.constrained == "partial") {
      hX.misreport.sensitive <- plogis(cbind(x.control, 1) %*% par.control, log.p = TRUE)
      hX.truthful.sensitive <- plogis(cbind(x.control, 1) %*% par.control, log.p = TRUE)
      hX.truthful.nonsensitive <- plogis(cbind(x.control, 0) %*% par.control, log.p = TRUE)
    }

    if(control.constrained == "full") {
      hX.misreport.sensitive <- plogis(x.control %*% par.control, log.p = TRUE)
      hX.truthful.sensitive <- plogis(x.control %*% par.control, log.p = TRUE)
      hX.truthful.nonsensitive <- plogis(x.control %*% par.control, log.p = TRUE)
    }

    hX.misreport.sensitive <- dbinom((y - treat * as.numeric(sensitive.response == 1)), size = J, prob = exp(hX.misreport.sensitive), log = TRUE)
    hX.truthful.sensitive <- dbinom((y - treat * as.numeric(sensitive.response == 1)), size = J, prob = exp(hX.truthful.sensitive), log = TRUE)
    hX.truthful.nonsensitive <- dbinom((y - treat * as.numeric(sensitive.response == 0)), size = J, prob = exp(hX.truthful.nonsensitive), log = TRUE)

    
    if(outcome.model != "none") {
      if(outcome.constrained == TRUE) {
        if(outcome.model %in% c("logistic", "binomial", "betabinomial")) {
          fX.misreport.sensitive <- plogis(cbind(x.outcome, 1) %*% par.outcome, log.p = TRUE)
          fX.truthful.sensitive <- plogis(cbind(x.outcome, 1) %*% par.outcome, log.p = TRUE)
          fX.truthful.nonsensitive <- plogis(cbind(x.outcome, 0) %*% par.outcome, log.p = TRUE)
        }
      } else {
        if(outcome.model %in% c("logistic", "binomial", "betabinomial")) {
          fX.misreport.sensitive <- plogis(cbind(x.outcome, 1, 1) %*% par.outcome, log.p = TRUE)
          fX.truthful.sensitive <- plogis(cbind(x.outcome, 0, 1) %*% par.outcome, log.p = TRUE)
          fX.truthful.nonsensitive <- plogis(cbind(x.outcome, 0, 0) %*% par.outcome, log.p = TRUE)
        }
      }
    } else {
        fX.misreport.sensitive <- rep(0, length(y))
        fX.truthful.sensitive <- rep(0, length(y))
        fX.truthful.nonsensitive <- rep(0, length(y))
    }

    if(outcome.model == "logistic") {
      fX.misreport.sensitive <- dbinom(o, size = 1, prob = exp(fX.misreport.sensitive), log = TRUE)
      fX.truthful.sensitive <- dbinom(o, size = 1, prob = exp(fX.truthful.sensitive), log = TRUE)
      fX.truthful.nonsensitive <- dbinom(o, size = 1, prob = exp(fX.truthful.nonsensitive), log = TRUE)
    } else if(outcome.model == "binomial") {
      fX.misreport.sensitive <- dbinom(o, size = trials, prob = exp(fX.misreport.sensitive), log = TRUE)
      fX.truthful.sensitive <- dbinom(o, size = trials, prob = exp(fX.truthful.sensitive), log = TRUE)
      fX.truthful.nonsensitive <- dbinom(o, size = trials, prob = exp(fX.truthful.nonsensitive), log = TRUE)
    } else if(outcome.model == "betabinomial") {
      fX.misreport.sensitive <- VGAM::dbetabinom(o, size = trials, prob = exp(fX.misreport.sensitive), rho = par.outcome.aux["rho"], log = TRUE)
      fX.truthful.sensitive <- VGAM::dbetabinom(o, size = trials, prob = exp(fX.truthful.sensitive), rho = par.outcome.aux["rho"], log = TRUE)
      fX.truthful.nonsensitive <- VGAM::dbetabinom(o, size = trials, prob = exp(fX.truthful.nonsensitive), rho = par.outcome.aux["rho"], log = TRUE)
    }

    
    if(sensitive.response == 1) {
      gX.misreport.sensitive <- plogis(x.sensitive %*% par.sensitive, log.p = TRUE)
      gX.truthful.sensitive <- plogis(x.sensitive %*% par.sensitive, log.p = TRUE)
      gX.truthful.nonsensitive <- log1p(-exp(plogis(x.sensitive %*% par.sensitive, log.p = TRUE))) 
    } else {
      gX.misreport.sensitive <- log1p(-exp(plogis(x.sensitive %*% par.sensitive, log.p = TRUE))) 
      gX.truthful.sensitive <- log1p(-exp(plogis(x.sensitive %*% par.sensitive, log.p = TRUE))) 
      gX.truthful.nonsensitive <- plogis(x.sensitive %*% par.sensitive, log.p = TRUE)
    }

    
    if(include.misreport.treatment == TRUE) {
      lX.misreport.sensitive <- plogis(cbind(x.misreport, treat) %*% par.misreport, log.p = TRUE)
      lX.truthful.sensitive <- log1p(-exp(plogis(cbind(x.misreport, treat) %*% par.misreport, log.p = TRUE))) 
      lX.truthful.nonsensitive <- log(rep(1, length(y))) 
    } else {
      lX.misreport.sensitive <- plogis(x.misreport %*% par.misreport, log.p = TRUE)
      lX.truthful.sensitive <- log1p(-exp(plogis(x.misreport %*% par.misreport, log.p = TRUE))) 
      lX.truthful.nonsensitive <- log(rep(1, length(y))) 
    }

    w[, 1] <- lX.misreport.sensitive + gX.misreport.sensitive + hX.misreport.sensitive + fX.misreport.sensitive
    w[, 2] <- lX.truthful.sensitive +     gX.truthful.sensitive +     hX.truthful.sensitive +     fX.truthful.sensitive
    w[, 3] <- lX.truthful.nonsensitive +  gX.truthful.nonsensitive +  hX.truthful.nonsensitive +  fX.truthful.nonsensitive

    w[respondentType == "Misreport sensitive", 1] <- log(1)
    w[respondentType == "Misreport sensitive", 2] <- log(0)
    w[respondentType == "Misreport sensitive", 3] <- log(0)

    w[respondentType == "Truthful sensitive", 1] <- log(0)
    w[respondentType == "Truthful sensitive", 2] <- log(1)
    w[respondentType == "Truthful sensitive", 3] <- log(0)

    w[respondentType == "Non-sensitive", 1] <- log(0)
    w[respondentType == "Non-sensitive", 2] <- log(0)
    w[respondentType == "Non-sensitive", 3] <- log(1)

    w[respondentType == "Non-sensitive or misreport sensitive", 2] <- log(0)

    denominator <- apply(w, 1, function(x) logAdd(logAdd(x[1], x[2]), x[3]))

    w[, 1] <- w[, 1] - denominator
    w[, 2] <- w[, 2] - denominator
    w[, 3] <- w[, 3] - denominator

    w[respondentType == "Misreport sensitive", 1] <- log(1)
    w[respondentType == "Misreport sensitive", 2] <- log(0)
    w[respondentType == "Misreport sensitive", 3] <- log(0)

    w[respondentType == "Truthful sensitive", 1] <- log(0)
    w[respondentType == "Truthful sensitive", 2] <- log(1)
    w[respondentType == "Truthful sensitive", 3] <- log(0)

    w[respondentType == "Non-sensitive", 1] <- log(0)
    w[respondentType == "Non-sensitive", 2] <- log(0)
    w[respondentType == "Non-sensitive", 3] <- log(1)

    w[respondentType == "Non-sensitive or misreport sensitive", 2] <- log(0)

    

    
    log.lik[respondentType == "Non-sensitive or misreport sensitive"] <- apply(data.frame(lX.truthful.nonsensitive[respondentType == "Non-sensitive or misreport sensitive"] +
                                                                                             gX.truthful.nonsensitive[respondentType == "Non-sensitive or misreport sensitive"] +
                                                                                             hX.truthful.nonsensitive[respondentType == "Non-sensitive or misreport sensitive"] +
                                                                                             fX.truthful.nonsensitive[respondentType == "Non-sensitive or misreport sensitive"],
                                                                                             lX.misreport.sensitive[respondentType == "Non-sensitive or misreport sensitive"] +
                                                                                             gX.misreport.sensitive[respondentType == "Non-sensitive or misreport sensitive"] +
                                                                                             hX.misreport.sensitive[respondentType == "Non-sensitive or misreport sensitive"] +
                                                                                             fX.misreport.sensitive[respondentType == "Non-sensitive or misreport sensitive"]),
                                                                                  1, function(x) logAdd(x[1], x[2]))

    
    log.lik[respondentType == "Truthful sensitive"] <- lX.truthful.sensitive[respondentType == "Truthful sensitive"] +
                                                       gX.truthful.sensitive[respondentType == "Truthful sensitive"] +
                                                       hX.truthful.sensitive[respondentType == "Truthful sensitive"] +
                                                       fX.truthful.sensitive[respondentType == "Truthful sensitive"]

    
    log.lik[respondentType == "Non-sensitive"] <- lX.truthful.nonsensitive[respondentType == "Non-sensitive"] +
                                                  gX.truthful.nonsensitive[respondentType == "Non-sensitive"] +
                                                  hX.truthful.nonsensitive[respondentType == "Non-sensitive"] +
                                                  fX.truthful.nonsensitive[respondentType == "Non-sensitive"]

    
    log.lik[respondentType == "Misreport sensitive"] <- lX.misreport.sensitive[respondentType == "Misreport sensitive"] +
                                                           gX.misreport.sensitive[respondentType == "Misreport sensitive"] +
                                                           hX.misreport.sensitive[respondentType == "Misreport sensitive"] +
                                                           fX.misreport.sensitive[respondentType == "Misreport sensitive"]

  }

  if(model.misreport == FALSE) {

    
    if(control.constrained == "none") {
      hX.1 <- plogis(cbind(x.control, 1) %*% par.control)
      hX.0 <- plogis(cbind(x.control, 0) %*% par.control)
    }

    if(control.constrained == "full") {
      hX.1 <- plogis(x.control %*% par.control)
      hX.0 <- plogis(x.control %*% par.control)
    }

    hX.1 <- dbinom((y - treat), size = J, prob = hX.1, log = TRUE)
    hX.0 <- dbinom(y, size = J, prob = hX.0, log = TRUE)

    
    if(outcome.model %in% c("logistic", "binomial", "betabinomial")) {
      fX.1 <- plogis(cbind(x.outcome, 1) %*% par.outcome)
      fX.0 <- plogis(cbind(x.outcome, 0) %*% par.outcome)
    } else {
      fX.1 <- rep(0, length(y))
      fX.0 <- rep(0, length(y))
    }

    if(outcome.model == "logistic") {
      fX.1 <- dbinom(o, size = 1, prob = fX.1, log = TRUE)
      fX.0 <- dbinom(o, size = 1, prob = fX.0, log = TRUE)
    } else if(outcome.model == "binomial") {
      fX.1 <- dbinom(o, size = trials, prob = fX.1, log = TRUE)
      fX.0 <- dbinom(o, size = trials, prob = fX.0, log = TRUE)
    } else if(outcome.model == "betabinomial") {
      fX.1 <- VGAM::dbetabinom(o, size = trials, prob = fX.1, rho = par.outcome.aux["rho"], log = TRUE)
      fX.0 <- VGAM::dbetabinom(o, size = trials, prob = fX.0, rho = par.outcome.aux["rho"], log = TRUE)
    }
    
    gX.1 <- plogis(x.sensitive %*% par.sensitive, log.p = TRUE)
    gX.0 <- log(1 - exp(gX.1))

    w[, 1] <- gX.1 + hX.1 + fX.1
    w[, 2] <- gX.0 + hX.0 + fX.0

    w[respondentType == "1", 1] <- log(1)
    w[respondentType == "1", 2] <- log(0)

    w[respondentType == "0", 1] <- log(0)
    w[respondentType == "0", 2] <- log(1)

    denominator <- apply(w, 1, function(x) logAdd(x[1], x[2]))

    w[, 1] <- w[, 1] - denominator
    w[, 2] <- w[, 2] - denominator

    w[respondentType == "1", 1] <- log(1)
    w[respondentType == "1", 2] <- log(0)

    w[respondentType == "0", 1] <- log(0)
    w[respondentType == "0", 2] <- log(1)

    log.lik[respondentType == "0"] <- gX.0[respondentType == "0"] +
                                      hX.0[respondentType == "0"] +
                                      fX.0[respondentType == "0"]

    log.lik[respondentType == "1"] <- gX.1[respondentType == "1"] +
                                      hX.1[respondentType == "1"] +
                                      fX.1[respondentType == "1"]

    log.lik[respondentType == "0 or 1"] <- apply(data.frame(gX.1[respondentType == "0 or 1"] +
                                                            hX.1[respondentType == "0 or 1"] +
                                                            fX.1[respondentType == "0 or 1"],
                                                            gX.0[respondentType == "0 or 1"] +
                                                            hX.0[respondentType == "0 or 1"] +
                                                            fX.0[respondentType == "0 or 1"]),
                                                 1, function(x) logAdd(x[1], x[2]))

    log.lik[respondentType == "0 or 1"] <- log(exp(gX.1[respondentType == "0 or 1"] +
                                                   hX.1[respondentType == "0 or 1"] +
                                                   fX.1[respondentType == "0 or 1"]) +
                                               exp(gX.0[respondentType == "0 or 1"] +
                                                   hX.0[respondentType == "0 or 1"] +
                                                   fX.0[respondentType == "0 or 1"]))
  }

  return(list(w = w, ll = sum(weight * log.lik)))
}

listExperiment <- function(formula, data, treatment, J, 
                           direct = NULL, sensitive.response = NULL,
                           outcome = NULL, outcome.trials = NULL,
                           outcome.model = "logistic",
                           outcome.constrained = TRUE,
                           control.constrained = "none",
                           include.misreport.treatment = TRUE,
                           weights = NULL, se = TRUE, tolerance = 1E-8,
                           max.iter = 10000, n.runs = 1, verbose = TRUE,
                           get.data = FALSE,
                           par.control = NULL, par.sensitive = NULL,
                           par.misreport = NULL, par.outcome = NULL,
                           par.outcome.aux = NULL,
                           formula.control = NULL, formula.sensitive = NULL,
                           formula.misreport = NULL, formula.outcome = NULL,
                           get.boot = 0, ...) {

  function.call <- match.call(expand.dots = FALSE)

  if(missing(data)) data <- environment(formula)
  mf <- match.call(expand.dots = FALSE)

  m <- match(c("formula", "data", "na.action"), names(mf), 0L)
  mf <- mf[c(1L, m)]
  mf$drop.unused.levels <- TRUE
  mf$na.action <- "na.pass"
  mf[[1]] <- quote(model.frame)
  mt <- attr(eval(mf, parent.frame()), "terms")
  xlevels.formula <- .getXlevels(attr(eval(mf, parent.frame()), "terms"), eval(mf, parent.frame()))

  if(!is.null(formula.control)) {
    mf.control <- mf
    mf.control$formula <- formula.control
    xlevels.formula.control <- .getXlevels(attr(eval(mf.control, parent.frame()), "terms"), eval(mf.control, parent.frame()))
    mf.control <- eval(mf.control, parent.frame())
    x.control <- model.matrix(attr(mf.control, "terms"), data = mf.control)
  } else {
    formula.control <- as.formula(mf$formula)
    xlevels.formula.control <- xlevels.formula
    mf.control <- eval(mf, parent.frame())
    x.control <- model.matrix(attr(mf.control, "terms"), data = mf.control)
  }

  if(!is.null(formula.sensitive)) {
    mf.sensitive <- mf
    mf.sensitive$formula <- formula.sensitive
    xlevels.formula.sensitive <- .getXlevels(attr(eval(mf.sensitive, parent.frame()), "terms"), eval(mf.sensitive, parent.frame()))
    mf.sensitive <- eval(mf.sensitive, parent.frame())
    x.sensitive <- model.matrix(attr(mf.sensitive, "terms"), data = mf.sensitive)
  } else {
    formula.sensitive <- as.formula(mf$formula)
    xlevels.formula.sensitive <- xlevels.formula
    mf.sensitive <- eval(mf, parent.frame())
    x.sensitive <- model.matrix(attr(mf.sensitive, "terms"), data = mf.sensitive)
  }

  if(!is.null(formula.misreport)) {
    mf.misreport <- mf
    mf.misreport$formula <- formula.misreport
    xlevels.formula.misreport <- .getXlevels(attr(eval(mf.misreport, parent.frame()), "terms"), eval(mf.misreport, parent.frame()))
    mf.misreport <- eval(mf.misreport, parent.frame())
    x.misreport <- model.matrix(attr(mf.misreport, "terms"), data = mf.misreport)
  } else {
    formula.misreport <- as.formula(mf$formula)
    xlevels.formula.misreport <- xlevels.formula
    mf.misreport <- eval(mf, parent.frame())
    x.misreport <- model.matrix(attr(mf.misreport, "terms"), data = mf.misreport)
  }

  if(!is.null(formula.outcome)) {
    mf.outcome <- mf
    mf.outcome$formula <- formula.outcome
    xlevels.formula.outcome <- .getXlevels(attr(eval(mf.outcome, parent.frame()), "terms"), eval(mf.outcome, parent.frame()))
    mf.outcome <- eval(mf.outcome, parent.frame())
    x.outcome <- model.matrix(attr(mf.outcome, "terms"), data = mf.outcome)
   } else {
    formula.outcome <- as.formula(mf$formula)
    xlevels.formula.outcome <- xlevels.formula
    mf.outcome <- eval(mf, parent.frame())
    x.outcome <- model.matrix(attr(mf.outcome, "terms"), data = mf.outcome)
  }

  mf <- eval(mf, parent.frame())
  y <- model.response(mf, type = "any")
  treat <- data[, paste(treatment)]

  xlevels <- c(xlevels.formula,
               xlevels.formula.control,
               xlevels.formula.sensitive,
               xlevels.formula.misreport,
               xlevels.formula.outcome)
  xlevels <- xlevels[-which(duplicated(xlevels))]


  
  x.control.na <- apply(x.control, 1, function(X) all(!is.na(X)))
  x.sensitive.na <- apply(x.sensitive, 1, function(X) all(!is.na(X)))
  x.misreport.na <- apply(x.misreport, 1, function(X) all(!is.na(X)))
  x.outcome.na <- apply(x.outcome, 1, function(X) all(!is.na(X)))

  y.na <- !is.na(y)
  treat.na <- !is.na(treat)

  if(!is.null(direct)) {
    d <- data[, paste(direct)]
    d.na <- !is.na(d)
    model.misreport <- TRUE
  } else {
    model.misreport <- FALSE
    d <- rep(NA, length(y))
    d.na <- rep(TRUE, length(y))
  }

  if(!is.null(outcome) & outcome.model %in% c("logistic")) {
    o <- data[, paste(outcome)]
    trials <- rep(NA, length(y))
    o.na <- !is.na(o)
  } else if(!is.null(outcome) & outcome.model %in% c("binomial", "betabinomial")) {
    o <- data[, paste(outcome)]
    trials <- data[, paste(outcome.trials)]
    o.na <- !is.na(o) & !is.na(trials)
  } else {
    o <- rep(NA, length(y))
    trials <- rep(NA, length(y))
    o.na <- rep(TRUE, length(y))
    outcome.model <- "none"
  }

  if(!is.null(weights)) {
    weight <- data[, paste(weights)]
    weight.na <- !is.na(weight)
  } else {
    weight <- rep(1, length(y))
    weight.na <- !is.na(weight)
  }

  y <- y[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na]
  x.control <- as.matrix(x.control[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na, , drop = FALSE])
  x.sensitive <- as.matrix(x.sensitive[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na, , drop = FALSE])
  x.outcome <- as.matrix(x.outcome[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na, , drop = FALSE])
  x.misreport <- as.matrix(x.misreport[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na, , drop = FALSE])
  treat <- treat[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na]
  d <- d[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na]
  o <- o[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na]
  trials <- trials[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na]
  weight <- weight[y.na & x.control.na & x.sensitive.na & x.outcome.na & x.misreport.na & treat.na & d.na & weight.na]

  n <- nrow(x.control)

  if(get.boot > 0) {
    boot.sample <- sample(1:length(weight), prob = weight, replace = TRUE)
    y <- as.matrix(y)[boot.sample, , drop = FALSE]
    x.control <- as.matrix(x.control)[boot.sample, , drop = FALSE]
    x.sensitive <- as.matrix(x.sensitive)[boot.sample, , drop = FALSE]
    x.outcome <- as.matrix(x.outcome)[boot.sample, , drop = FALSE]
    x.misreport <- as.matrix(x.misreport)[boot.sample, , drop = FALSE]
    treat <- as.matrix(treat)[boot.sample]
    d <- as.matrix(d)[boot.sample, , drop = FALSE]
    o <- as.matrix(o)[boot.sample, , drop = FALSE]
    trials <- as.matrix(trials)[boot.sample, , drop = FALSE]
    weight <- rep(1, length(y))
    se <- FALSE
  }

  respondentType <- rep(as.character(NA), length(y))

  if(model.misreport == TRUE) {

    
    respondentType[treat == 0 & d != sensitive.response] <- "Non-sensitive or misreport sensitive"

    
    respondentType[treat == 0 & d == sensitive.response] <- "Truthful sensitive"

    
    if(sensitive.response == 1) respondentType[treat == 1 & y == 0 & d != sensitive.response] <- "Non-sensitive"
    if(sensitive.response == 0) respondentType[treat == 1 & y == (J + 1) & d != sensitive.response] <- "Non-sensitive"

    
    respondentType[treat == 1 & y > 0 & y < (J + 1) & d != sensitive.response] <- "Non-sensitive or misreport sensitive"

    
    if(sensitive.response == 1) respondentType[treat == 1 & y == (J + 1) & d != sensitive.response] <- "Misreport sensitive"
    if(sensitive.response == 0) respondentType[treat == 1 & y == 0 & d != sensitive.response] <- "Misreport sensitive"

    
    if(sensitive.response == 1) respondentType[treat == 1 & y == (J + 1) & d == sensitive.response] <- "Truthful sensitive"
    if(sensitive.response == 0) respondentType[treat == 1 & y == 0 & d == sensitive.response] <- "Truthful sensitive"

    
    if(sensitive.response == 1) respondentType[treat == 1 & y == 0 & d == sensitive.response] <- "Violates assumption"
    if(sensitive.response == 0) respondentType[treat == 1 & y == (J + 1) & d == sensitive.response] <- "Violates assumption"

    
    respondentType[treat == 1 & y > 0 & y < (J + 1) & d == sensitive.response] <- "Truthful sensitive"
  } else {
    
    respondentType[treat == 1 & y > 0 & y < (J + 1)] <- "0 or 1"

    
    respondentType[treat == 0] <- "0 or 1"

    
    respondentType[treat == 1 & y == 0] <- "0"

    
    respondentType[(treat == 1 & y == (J + 1))] <- "1"
  }

  if("Violates assumption" %in% respondentType) {
    stop("\nSome observations violate the monotonicity assumption.")
  }

  if(model.misreport == TRUE) {
    w <- as.matrix(data.frame(as.numeric(respondentType %in% c("Non-sensitive or misreport sensitive", "Misreport sensitive")),
                              as.numeric(respondentType == "Truthful sensitive"),
                              as.numeric(respondentType %in% c("Non-sensitive or misreport sensitive", "Non-sensitive"))))
    w <- w / apply(w, 1, sum)
    colnames(w) <- c("Misreport sensitive", "Truthful sensitive", "Non-sensitive")
  } else {
    w <- as.matrix(data.frame(as.numeric(respondentType %in% c("1", "0 or 1")),
                              as.numeric(respondentType %in% c("0", "0 or 1"))))
    w <- w / apply(w, 1, sum)
    colnames(w) <- c("1", "0")
  }

  w <- log(w)

  if(get.data == TRUE) {

    estep.out <- estep(y = y, w = w, x.control = x.control, x.sensitive = x.sensitive, x.outcome = x.outcome, x.misreport = x.misreport, treat = treat, J = J,
                       par.sensitive = par.sensitive, par.control = par.control, par.outcome = par.outcome, par.outcome.aux = par.outcome.aux, par.misreport = par.misreport,
                       d = d, sensitive.response = sensitive.response, model.misreport = model.misreport,
                       o = o, trials = trials, outcome.model = outcome.model, 
                       weight = weight, respondentType = respondentType,
                       outcome.constrained = outcome.constrained,
                       control.constrained = control.constrained,
                       include.misreport.treatment = include.misreport.treatment)

    return(list(w = estep.out$w,
                ll = estep.out$ll,
                x.control = x.control,
                x.sensitive = x.sensitive,
                x.misreport = x.misreport,
                x.outcome = x.outcome))

  }

  if(model.misreport == TRUE) {
    if(include.misreport.treatment == TRUE) par.misreport <- rep(0, ncol(x.misreport) + 1) 
    if(include.misreport.treatment == FALSE) par.misreport <- rep(0, ncol(x.misreport))
  } else {
    par.misreport <- NULL
  }

  if(is.null(par.sensitive)) par.sensitive <- rep(0, ncol(x.sensitive))

  if(is.null(par.control)) {
    if(control.constrained == "none" & model.misreport == FALSE) {
      par.control <- rep(0, ncol(x.control) + 1)
    } else if(control.constrained == "none" & model.misreport == TRUE) {
      par.control <- rep(0, ncol(x.control) + 2)
    } else if(control.constrained == "partial" & model.misreport == FALSE) {
      stop("If not modeling misreporting, set argument control.constrained to 'none' or 'full'")
    } else if(control.constrained == "partial" & model.misreport == TRUE) {
      par.control <- rep(0, ncol(x.control) + 1)
    } else if(control.constrained == "full") {
      par.control <- rep(0, ncol(x.control))
    }
  }

  if(is.null(par.outcome)) {
    if(outcome.model != "none") {
      if(outcome.constrained == TRUE) par.outcome <- rep(0, ncol(x.outcome) + 1)
      if(outcome.constrained == FALSE) par.outcome <- rep(0, ncol(x.outcome) + 2)
    } else {
      par.outcome <- NULL
    }
  }

  if(is.null(par.outcome.aux)) {
    if(outcome.model %in% c("none", "logistic")) {
      par.outcome.aux <- NULL
    } else if(outcome.model == "betabinomial") {
      par.outcome.aux <- list(rho = 0)
    }
  }
  
  runs <- list()

  
  for(j in 1:n.runs) {

    logLikelihood <- rep(as.numeric(NA), max.iter)

    while(TRUE) {
      par.control <- rnorm(length(par.control), 0.25)
      par.sensitive <- rnorm(length(par.sensitive), 0.25)
      if(model.misreport == TRUE) par.misreport <- rnorm(length(par.misreport), 0.25)
      if(outcome.model != "none") par.outcome <- rnorm(length(par.outcome), 0.25)
      if(outcome.model != "none" & length(par.outcome.aux) > 0) par.outcome.aux <- abs(rnorm(length(par.outcome.aux), 0.25))

      templl <- estep(y = y, w = w, x.control = x.control, x.sensitive = x.sensitive, x.outcome = x.outcome, x.misreport = x.misreport, treat = treat, J = J,
                      par.sensitive = par.sensitive, par.control = par.control, par.outcome = par.outcome, par.outcome.aux = par.outcome.aux, par.misreport = par.misreport,
                      d = d, sensitive.response = sensitive.response, model.misreport = model.misreport,
                      o = o, trials = trials, outcome.model = outcome.model, 
                      weight = weight, respondentType = respondentType,
                      outcome.constrained = outcome.constrained,
                      control.constrained = control.constrained,
                      include.misreport.treatment = include.misreport.treatment)$ll
      templl

      if(!is.nan(templl) & templl > -Inf) break()
    }

    for(i in 1:max.iter) { 

      estep.out <- estep(y = y, w = w, x.control = x.control, x.sensitive = x.sensitive, x.outcome = x.outcome, x.misreport = x.misreport, treat = treat, J = J,
                         par.sensitive = par.sensitive, par.control = par.control, par.outcome = par.outcome, par.outcome.aux = par.outcome.aux, par.misreport = par.misreport,
                         d = d, sensitive.response = sensitive.response, model.misreport = model.misreport,
                         o = o, trials = trials, outcome.model = outcome.model, 
                         weight = weight, respondentType = respondentType,
                         outcome.constrained = outcome.constrained,
                         control.constrained = control.constrained,
                         include.misreport.treatment = include.misreport.treatment)

      w <- estep.out$w
      logLikelihood[i] <- estep.out$ll

      if(i > 1 & verbose == TRUE & get.boot == 0) {
        cat("\r\rRun:", paste0(j, "/", n.runs), "Iter:", i,
            "llik:", sprintf("%.2f", logLikelihood[i]),
            "llik change:", sprintf("%.8f", (logLikelihood[i] - logLikelihood[i-1])),
            "(tol =", paste0(as.character(tolerance), ")       "))
      }

      if(i > 1 & verbose == TRUE & get.boot > 0) {
        cat("\r\rBoot:", get.boot, "Run:", paste0(j, "/", n.runs), "Iter:", i,
            "llik:", sprintf("%.2f", logLikelihood[i]),
            "llik change:", sprintf("%.8f", (logLikelihood[i] - logLikelihood[i-1])),
            "(tol =", paste0(as.character(tolerance), ")       "))
      }

      if(i > 1 && (logLikelihood[i] - logLikelihood[i - 1]) < 0) {
        stop("Log-likelihood increasing.")
      }
      if(i > 1 && (logLikelihood[i] - logLikelihood[i - 1]) < tolerance) {
        break()
      }

      par.sensitive <- mstepSensitive(y = y, treat = treat, x.sensitive = x.sensitive, w = w,
                                      d = d, sensitive.response = sensitive.response,
                                      weight = weight, model.misreport = model.misreport)

      par.control <- mstepControl(y = y, J = J, treat = treat, x.control = x.control, w = w,
                                  d = d, sensitive.response = sensitive.response,
                                  weight = weight, model.misreport = model.misreport,
                                  control.constrained = control.constrained)

      if(outcome.model != "none") {
        outcome <- mstepOutcome(y = y, treat = treat, x.outcome = x.outcome, w = w,
                                d = d, sensitive.response = sensitive.response,
                                o = o, trials = trials, weight = weight,
                                model.misreport = model.misreport,
                                outcome.model = outcome.model,
                                outcome.constrained = outcome.constrained,
                                control.constrained = control.constrained)
        par.outcome <- outcome$coefs
        par.outcome.aux <- outcome$coefs.aux
      }

      if(model.misreport == TRUE) {
        par.misreport <- mstepMisreport(y = y, x.misreport = x.misreport,
                                              w = w, treat = treat,
                                              include.misreport.treatment = include.misreport.treatment,
                                              weight = weight)
      }
    }
  
    runs[[j]] <- list(logLikelihood = logLikelihood[i],
                      par.control = par.control,
                      par.sensitive = par.sensitive,
                      par.misreport = par.misreport,
                      par.outcome = par.outcome,
                      par.outcome.aux = par.outcome.aux)

  }

  if(verbose == TRUE) cat("\n")

  max.ll <- which(sapply(runs, function(X) X$logLikelihood) == max(sapply(runs, function(X) X$logLikelihood)))
  llik <- runs[[max.ll]]$logLikelihood
  par.control <- runs[[max.ll]]$par.control
  par.sensitive <- runs[[max.ll]]$par.sensitive
  par.misreport <- runs[[max.ll]]$par.misreport
  par.outcome <- runs[[max.ll]]$par.outcome
  par.outcome.aux <- runs[[max.ll]]$par.outcome.aux

  par <- c(par.control, par.sensitive, par.misreport, par.outcome, par.outcome.aux)
  num <- c(length(par.control), length(par.sensitive), length(par.misreport), length(par.outcome), length(par.outcome.aux))

  llik.wrapper <- function(par, num, y, w,
                           x.control, x.sensitive, x.outcome, x.misreport, treat, J,
                           d, sensitive.response, model.misreport,
                           o, trials, outcome.model,
                           weight, respondentType,
                           outcome.constrained,
                           control.constrained,
                           include.misreport.treatment) {
    par.control <- par[1:num[1]]
    par.sensitive <- par[(num[1]+1):sum(num[1:2])]
    if(model.misreport == TRUE) {
      par.misreport <- par[(sum(num[1:2])+1):sum(num[1:3])]
    } else{
      par.misreport <- NULL
    }
    if(outcome.model != "none") {
      par.outcome <- par[(sum(num[1:3])+1):sum(num[1:4])]
      if(outcome.model %in% c("betabinomial", "linear")) {
        par.outcome.aux <- par[(sum(num[1:4])+1):sum(num[1:5])]
      } else {
        par.outcome.aux <- NULL
      }
    } else {
      par.outcome <- NULL
    }

    llik <- estep(y = y, w = w, x.control = x.control, x.sensitive = x.sensitive, x.outcome = x.outcome, x.misreport = x.misreport, treat = treat, J = J,
                  par.sensitive = par.sensitive, par.control = par.control, par.outcome = par.outcome, par.outcome.aux = par.outcome.aux, par.misreport = par.misreport,
                  d = d, sensitive.response = sensitive.response, model.misreport = model.misreport,
                  o = o, trials = trials, outcome.model = outcome.model,
                  weight = weight, respondentType = respondentType,
                  outcome.constrained = outcome.constrained,
                  control.constrained = control.constrained,
                  include.misreport.treatment)$ll
    return(llik)
  }

  if(se == TRUE & all(weight == 1)) {

    num <- c(length(par.control),
             length(par.sensitive),
             length(par.misreport),
             length(par.outcome),
             length(par.outcome.aux))

    hess <- numDeriv::hessian(llik.wrapper, c(par.control, par.sensitive, par.misreport, par.outcome, par.outcome.aux),
                              num = num, J = J, y = y, w = w, treat = treat,
                              x.control = x.control, x.sensitive = x.sensitive, x.outcome = x.outcome, x.misreport = x.misreport,
                              d = d, sensitive.response = sensitive.response, model.misreport = model.misreport,
                              o = o, outcome.model = outcome.model,
                              outcome.constrained = outcome.constrained,
                              weight = weight,
                              respondentType = respondentType,
                              control.constrained = control.constrained,
                              include.misreport.treatment = include.misreport.treatment,
                              method.args = list(zero.tol = 1e-10))

    vcov.mle <- solve(-hess)
    se.mle <- sqrt(diag(vcov.mle))
    
    se.control <- se.mle[1:num[1]]
    names(se.control) <- names(par.control)
    
    se.sensitive <- se.mle[(num[1]+1):sum(num[1:2])]
    names(se.sensitive) <- names(par.sensitive)
    
    if(model.misreport == TRUE) {
      se.misreport <- se.mle[(sum(num[1:2])+1):sum(num[1:3])]
      names(se.misreport) <- names(par.misreport)
    } else {
      se.misreport <- NULL
    }
    if(outcome.model != "none") {
      se.outcome <- se.mle[(sum(num[1:3])+1):sum(num[1:4])]
      names(se.outcome) <- names(par.outcome)
      if(outcome.model %in% c("linear", "betabinomial")) {
        se.outcome.aux <- se.mle[(sum(num[1:4])+1):sum(num[1:5])]
        names(se.outcome.aux) <- names(par.outcome.aux)
      } else {
        se.outcome.aux <- NULL
      }
    } else {
     se.outcome <- NULL
     se.outcome.aux <- NULL
    }

  } else {
    se.control <- se.sensitive <- se.misreport <- se.outcome <- se.outcome.aux <- vcov.mle <- NULL
    if(se == TRUE) {
      warning("Standard errors are not implemented for models with survey weights.")
      se <- FALSE
    }
  }

  return.object <- list("par.control" = par.control,
                        "par.sensitive" = par.sensitive,
                        "par.misreport" = par.misreport,
                        "par.outcome" = par.outcome,
                        "par.outcome.aux" = par.outcome.aux,
                        "df" = n - length(c(par.control, par.sensitive, par.misreport, par.outcome, par.outcome.aux)),
                        "se.sensitive" = se.sensitive,
                        "se.control" = se.control,
                        "se.misreport" = se.misreport,
                        "se.outcome" = se.outcome,
                        "se.outcome.aux" = se.outcome.aux,
                        "vcov.mle" = vcov.mle,
                        "w" = w,
                        "data" = data,
                        "direct" = direct,
                        "treatment" = treatment,
                        "model.misreport" = model.misreport,
                        "outcome.model" = outcome.model,
                        "outcome.constrained" = outcome.constrained,
                        "control.constrained" = control.constrained,
                        "include.misreport.treatment" = include.misreport.treatment,
                        "weights" = weights,
                        "formula" = formula,
                        "formula.control" = formula.control,
                        "formula.sensitive" = formula.sensitive,
                        "formula.misreport" = formula.misreport,
                        "formula.outcome" = formula.outcome,
                        "sensitive.response" = sensitive.response,
                        "xlevels" = xlevels,
                        "llik" = llik,
                        "n" = n,
                        "J" = J,
                        "se" = se,
                        "runs" = runs,
                        "call" = function.call,
                        "boot" = FALSE)

  class(return.object) <- "listExperiment"
  return(return.object)

}


bootListExperiment <- function(formula, data, treatment, J, 
                               direct = NULL, sensitive.response = NULL,
                               outcome = NULL, outcome.trials = NULL, outcome.model = "logistic",
                               outcome.constrained = TRUE, control.constrained = "partial",
                               include.misreport.treatment = TRUE,
                               weights = NULL, se = TRUE, tolerance = 1E-8, max.iter = 5000,
                               n.runs = 10, verbose = TRUE, get.data = FALSE,
                               par.control = NULL, par.sensitive = NULL, par.misreport = NULL,
                               par.outcome = NULL, par.outcome.aux = NULL,
                               formula.control = NULL, formula.sensitive = NULL,
                               formula.misreport = NULL, formula.outcome = NULL,
                               boot.iter = 1000, parallel = FALSE, n.cores = 2, cluster = NULL) {

  function.call <- match.call()
  args.call <- as.list(function.call)[-1]
  args.call$se <- FALSE
  args.call$get.boot <- 1

  args.call <- lapply(args.call, eval)

  data <- args.call$data
  args.call$data <- as.name("data")
  
  if(parallel == FALSE) {
    boot.out <- list()
    for(i in 1:boot.iter) {   
      args.call$get.boot <- i
      boot.out[[i]] <- do.call(listExperiment, args.call)
    }
  }

  if(parallel == TRUE) {
      
    args.call$verbose <- FALSE

    cat("Running bootstrap in parallel on ", n.cores, " cores/threads (", parallel::detectCores(), " available)...\n", sep = ""); Sys.sleep(0.2)

    if(!is.null(cluster)) cl <- cluster
    if(is.null(cluster)) cl <- parallel::makeCluster(n.cores)
    parallel::clusterExport(cl,
                  list("args.call", "data", "listExperiment", "logAdd", "estep",
                       "mstepControl", "mstepSensitive", "mstepMisreport", "mstepOutcome"),
                  envir = environment())

    boot.out <- parallel::parLapply(cl, 1:boot.iter, function(x) do.call(listExperiment, args.call))
    parallel::stopCluster(cl)
  }

  getPars <- function(varName) {
    X <- do.call(rbind, sapply(boot.out, function(x) x[varName]))
    cov.var <- cov(X)
    par.var <- colMeans(X)
    se.var <- as.vector(as.matrix(sqrt(diag(cov.var))))
    names(se.var) <- row.names(cov.var)
    return(list(par = par.var, se = se.var))
  }

  par.control <- getPars("par.control")$par
  se.control <- getPars("par.control")$se
  
  par.sensitive <- getPars("par.sensitive")$par
  se.sensitive <- getPars("par.sensitive")$se
  
  if(!is.null(boot.out[[1]]$par.misreport)) {
      par.misreport <- getPars("par.misreport")$par
      se.misreport <- getPars("par.misreport")$se
  } else {
    par.misreport <- se.misreport <- NULL
  }

  if(!is.null(boot.out[[1]]$par.outcome)) {
    par.outcome <- getPars("par.outcome")$par
    se.outcome <- getPars("par.outcome")$se
  } else {
    par.outcome <- se.outcome <- NULL
  }

  if(!is.null(boot.out[[1]]$outcome.model.aux)) {
    par.outcome <- getPars("par.outcome.aux")$par
    se.outcome <- getPars("par.outcome.aux")$se
  } else {
    par.outcome.aux <- se.outcome.aux <- NULL
  }

  se <- TRUE

  args.call$get.boot <- 0
  args.call$get.data <- TRUE
  args.call$par.control <- par.control
  args.call$par.sensitive <- par.sensitive
  args.call$par.misreport <- par.misreport
  args.call$par.outcome <- par.outcome
  args.call$par.outcome.aux <- par.outcome.aux
  llik <- do.call(listExperiment, args.call)$ll
  w <- do.call(listExperiment, args.call)$w

  return.object <- boot.out[[1]] 
  return.object$par.control <- par.control
  return.object$par.sensitive <- par.sensitive
  return.object$par.misreport <- par.misreport
  return.object$par.outcome <- par.outcome
  return.object$par.outcome.aux <- par.outcome.aux
  return.object$df <- return.object$n - length(c(par.control, par.sensitive, par.misreport, par.outcome, par.outcome.aux))
  return.object$se.control <- se.control
  return.object$se.sensitive <- se.sensitive
  return.object$se.misreport <- se.misreport
  return.object$se.outcome <- se.outcome
  return.object$se.outcome.aux <- se.outcome.aux
  return.object$vcov.model <- NULL
  return.object$data <- data
  return.object$se <- TRUE
  
  return.object$w <- w
  return.object$llik <- llik
  return.object$call <- function.call
  return.object$boot.iter <- boot.iter
  return.object$boot.out <- boot.out
  return.object$boot <- TRUE

  class(return.object) <- "listExperiment"
  return(return.object)

}


summary.listExperiment <- function(object, digits = 3, show.p = TRUE, model = c("control", "sensitive", "misreport", "outcome"), intercept = "top", ...) {

  Control <- Sensitive <- Misreport <- Outcome <- NULL

  # CONTROL-ITEMS SUB-MODEL
  par.control <- se.control <- z.control <- p.control <- star.control <- NULL
  Control <- data.frame(est = round(object$par.control, digits))
  if(object$se == TRUE) {
    Control <- data.frame(Control, se = round(object$se.control, digits))
    if(show.p == TRUE) Control <- data.frame(Control, p = round(2 * pnorm(abs(object$par.control/object$se.control), lower.tail = FALSE), digits))
    Control <- data.frame(Control, star = cut((2 * pnorm(abs(object$par.control/object$se.control), lower.tail = FALSE)),
                                              c(0, 0.001, 0.01, 0.05, 1), label = c("***", "** ", "*  ", "   ")))
    names(Control)[names(Control) == "star"] <- ""
  }
  if(intercept == "bottom") Control <- Control[c(2:nrow(Control), 1), ]

  # SENSITIVE-ITEM SUB-MODEL
  par.sensitive <- se.sensitive <- z.sensitive <- p.sensitive <- star.sensitive <- NULL
  Sensitive <- data.frame(est = round(object$par.sensitive, digits))
  if(object$se == TRUE) {
    Sensitive <- data.frame(Sensitive, se = round(object$se.sensitive, digits))
    if(show.p == TRUE) Sensitive <- data.frame(Sensitive, p = round(2 * pnorm(abs(object$par.sensitive/object$se.sensitive), lower.tail = FALSE), digits))
    Sensitive <- data.frame(Sensitive, star = cut((2 * pnorm(abs(object$par.sensitive/object$se.sensitive), lower.tail = FALSE)),
                                              c(0, 0.001, 0.01, 0.05, 1), label = c("***", "** ", "*  ", "   ")))
    names(Sensitive)[names(Sensitive) == "star"] <- ""
  }
  if(intercept == "bottom") Sensitive <- Sensitive[c(2:nrow(Sensitive), 1), ]

  # MISREPORT SUB-MODEL
  if(object$model.misreport == TRUE) {

    par.misreport <- se.misreport <- z.misreport <- p.misreport <- star.misreport <- NULL
    Misreport <- data.frame(est = round(object$par.misreport, digits))
    if(object$se == TRUE) {
      Misreport <- data.frame(Misreport, se = round(object$se.misreport, digits))
      if(show.p == TRUE) Misreport <- data.frame(Misreport, p = round(2 * pnorm(abs(object$par.misreport/object$se.misreport), lower.tail = FALSE), digits))
      Misreport <- data.frame(Misreport, star = cut((2 * pnorm(abs(object$par.misreport/object$se.misreport), lower.tail = FALSE)),
                                                c(0, 0.001, 0.01, 0.05, 1), label = c("***", "** ", "*  ", "   ")))
      names(Misreport)[names(Misreport) == "star"] <- ""
    }
    if(intercept == "bottom") Misreport <- Misreport[c(2:nrow(Misreport), 1), ]

  }

  # OUTCOME SUB-MODEL
  if(object$outcome.model != "none") {

    par.outcome <- se.outcome <- z.outcome <- p.outcome <- star.outcome <- NULL
    Outcome <- data.frame(est = round(object$par.outcome, digits))
    if(object$se == TRUE) {
      Outcome <- data.frame(Outcome, se = round(object$se.outcome, digits))
      if(show.p == TRUE) Outcome <- data.frame(Outcome, p = round(2 * pnorm(abs(object$par.outcome/object$se.outcome), lower.tail = FALSE), digits))
      Outcome <- data.frame(Outcome, star = cut((2 * pnorm(abs(object$par.outcome/object$se.outcome), lower.tail = FALSE)),
                                                c(0, 0.001, 0.01, 0.05, 1), label = c("***", "** ", "*  ", "   ")))
      names(Outcome)[names(Outcome) == "star"] <- ""
    }
    if(intercept == "bottom") Outcome <- Outcome[c(2:nrow(Outcome), 1), ]

  }

  cat("\nList experiment sub-models\n\n")

  cat("Call: ")
  print(object$call)

  if("control" %in% model) {
    cat("\nCONTROL ITEMS Pr(Y* = y)\n")
    print(Control); cat("---\n")
  }

  if("sensitive" %in% model) {
    cat("\nSENSITIVE ITEM Pr(Z* = 1)\n")
    print(Sensitive); cat("---\n")
  }

  if("misreport" %in% model & !is.null(Misreport)) {
    cat("\nMISREPORT Pr(U* = 1)\n")
    print(Misreport); cat("---\n")
  }

  if("outcome" %in% model & !is.null(Outcome)) {
    cat("\nOUTCOME\n")
    print(Outcome); cat("---\n")
  }

  if(object$se == TRUE) cat("* p < 0.05; ** p < 0.01; *** p < 0.001\n")

  if(object$boot == TRUE) {
    cat("\nStandard errors calculated by non-parametric bootstrap (", format(object$boot.iter, big.mark = ","), " draws).", sep = "")
  }

  cat("\nObservations:", format(object$n, big.mark = ","))
  # if(nrow(object$data)-object$n != 0)
  cat(" (", format(nrow(object$data)-object$n, big.mark = ","), " of ", format(nrow(object$data), big.mark = ","), " observations removed due to missingness)", sep = "")
  cat("\nLog-likelihood", object$llik)

}
